| Multivariate Bernstein Polynomials

> {-# LANGUAGE ViewPatterns #-}
> {-# OPTIONS_GHC -fno-warn-overlapping-patterns #-}

> module Bs (
>   Bs, ofPly, diff, 
>   lb, ub, lbubN, 
>   Scaled(Scaled), mkScaled, 
>   subdivl, subdivr, BsClass(..)) where

> import Aux (map', zipKeep, rapp, fact, Hand(..), (!:))
> import qualified Tree
> import Tree (Tree(..))
> import qualified Ply
> import Data.List (foldl1')
> import Numeric (showFFloat)

| The type of multivariate Bernstein polynomials with coefficients in
  @a@.

> newtype Bs a = B (Tree a) deriving (Eq,Show)


TODO make abstract

| A polynomial wraped together with a scale factor.

> data Scaled a = Scaled !(Bs a) !a deriving (Eq,Show)

> unB :: Bs t -> Tree t
> unB (B b) = b

> instance Num a => Num (Bs a) where
>   B (Lf x)       + B (Lf y)       = B $ Lf $! x+y
>   B x            + B (Br [])      = B x
>   B (Br [])      + B y            = B y
>   B (Br x)       + B (Br y)       = B $ Br $! zipKeep (|+|) x y
>   B x            + B (Br (y0:y')) = B $ Br (x|+|y0 !: y')
>   B (Br (x0:x')) + B y            = B $ Br (x0|+|y !: x')

>   negate (B a) = B $ fmap negate a

>   (*) = error "a"
>   abs = error "a"
>   signum = error "a"
>   fromInteger = error "a"

> instance Functor Bs where
>   fmap f (B p) = B $ fmap f p

> (|+|) :: Num a => Tree a -> Tree a -> Tree a
> p |+| q = let B r = B p + B q in r


* Integration

> integT :: Num a => Tree a -> Tree a -> Tree a
> integT c d@(Lf _) = Br $! [c, unB (B c + B d)]
> integT c (Br b) = Br $! integ'T c b

> integ :: Num a => Bs a -> Bs a -> Bs a
> integ c b = B (integT (unB c) (unB b))

> integ'T :: Num a => Tree a -> [Tree a] -> [Tree a]
> integ'T c [] = [c]
> integ'T c (b0:b') = c !: integ'T (unB (B c + B b0)) b'

> integ' :: Num t => t -> [t] -> [t]
> integ' c [] = [c]
> integ' c (b0:b') = c !: integ' (c + b0) b'

> (.*) :: Num a => a -> Bs a -> Bs a
> s .* p = fmap (s*) p


* Conversion from Power Basis

| Conversion from the Horner basis (in the first variable)

> ofPly1 :: Ply.Ply Double -> [Ply.Ply Double]
> ofPly1 = Ply.rec_integ [] (\c b r -> integ' (fmap ((fromInteger $ fact $ Ply.len b) *) c) r)

| Conversion from the Horner basis.

> ofPly :: Ply.Ply Double -> Bs Double
> ofPly (Ply.isConst -> Just c) = B $ Lf $! c
> ofPly b                       = B $ Br $! map' (unB . ofPly) (ofPly1 b)


* Bounding
 
| The scale factor, assuming the argument is levelled.

> scale :: Num a => Bs a -> Integer
> scale (B (Br (b0:b'))) = fact (length b') * scale (B b0)
> scale _ = 1

> mkScaled :: Fractional a => Bs a -> Scaled a
> mkScaled b = Scaled b (recip $ fromInteger $ scale b)

| A lower bound on the given polynomial.

> lb :: (Num a, Ord a) => Scaled a -> a
> lb (Scaled (B b) s) = Tree.fold 0 min b * s

| An upper bound on the given polynomial.

> ub :: (Num a, Ord a) => Scaled a -> a
> ub (Scaled (B b) s) = Tree.fold 0 max b * s

> corners :: Num a => Int -> Scaled a -> [a]
> corners k (Scaled (B b) s) = corners_ k b
>   where
>     corners_ 0     (Lf c) = [c*s]
>     corners_ (k+1) c@(Lf _) = let x = corners_ k c in x `seq` (x++x)
>     corners_ (k+1) c@(Br b) = corners_ k (head b) `rapp` corners_ k (last b)

| An upper bound on the lower bound.

> lbubN :: (Num a, Ord a) => Int -> [Scaled a] -> a
> lbubN n [] = error "lbubN.empty"
> lbubN n bs' = foldl1' min $ foldl1' (zipWith max) $ map' (corners n) bs'


* Differentiation and Differential Recursion

> diffT :: Num a => Tree a -> [Tree a] -> [Tree a]
> diffT _  []      = []
> diffT b0 (b1:b') = unB (B b1 - B b0) !: diffT b1 b' 

| Differentiation w.r.t the given dimension.

> diff :: Num a => Int -> Bs a -> Bs a
> diff 0     (B (Br (b0:b'))) = B $ Br $! diffT b0 b'
> diff (k+1) (B (Br b))       = B $ Br $! map' (unB . diff k . B) b
> diff _     _                = B $ Br []

> rec :: Num b => a -> (Bs b -> Bs b -> a -> a) -> Bs b -> a
> rec z _   (B (Br []))     = z
> rec z f p@(B (Lf c))      = let d = diff 0 p in f p d (rec z f d)
> rec z f p@(B (Br (b0:_))) = let d = diff 0 p in f (B b0) d (rec z f d)


* Subdivision

TODO wrap these in BsClass

| Subdivide into two (left part)

> subdivl :: Fractional a => Int -> Bs a -> Bs a
> subdivl 0     b = (0.5 ^ (Tree.width (unB b) - 1)) .* 
>                   (B $ rec (Br []) 
>                     (\c b r -> unB $ integ ((2 ^ Tree.width (unB b)) .* c) (B r)) b)
> subdivl (k+1) (B (Br b)) = B $ Br $! map' (unB . subdivl k . B) b
> subdivl (k+1) (B (Lf c)) = B $ Lf c

| Subdivide into two (right part)

> subdivr :: Fractional a => Int -> Bs a -> Bs a
> subdivr 0     (B b)      = B $ Tree.rev $ unB $ subdivl 0 $ B $ Tree.rev b
> subdivr (k+1) (B (Br b)) = B $ Br $! map' (unB . subdivr k . B) b
> subdivr (k+1) (B (Lf c)) = B $ Lf c

> subdivlS n (Scaled b s) = Scaled (subdivl n b) s
> subdivlS :: Fractional t => Int -> Scaled t -> Scaled t

> subdivrS :: Fractional t => Int -> Scaled t -> Scaled t
> subdivrS n (Scaled b s) = Scaled (subdivr n b) s


*** Sharpness

> le :: Ord a => Tree a -> Tree a -> Bool
> le (Lf x) (Lf y) = x <= y
> le (Br x) (Br y) = all id $ zipWith le x y
> le _      _      = error "Ouch! A branch fell on a leaf!"
>                    -- impossible with Bernstein polys obtained through ofPly

> sharpTreeL :: Ord a => Int -> Tree a -> Bool
> sharpTreeL (k+1) (Br b)       = all (sharpTreeL k) b
> sharpTreeL 0     (Br (b0:b')) = all (le b0) b'
> sharpTreeL _     _            = True

> sharpTreeR :: Ord a => Int -> Tree a -> Bool
> sharpTreeR (k+1) (Br b) = all (sharpTreeR k) b
> sharpTreeR 0     (Br b) = all (le $ last b) b
> sharpTreeR _      _     = True

> projTree :: Int -> Hand -> Tree a -> Tree a
> projTree _     _ (l@(Lf _)) = l
> projTree 0     L (Br b)     = head b
> projTree 0     R (Br b)     = last b
> projTree (k+1) h (Br b)     = Br $ map (projTree k h) b


TODO rename Bs to T and BsClass to Bs

> class BsClass b where
>   sharp :: Int -> b -> Maybe Hand
>   -- ^ Checks whether the polynomial is sharp w.r.t the given dimension.
>   proj :: Int -> Hand -> b -> b
>   -- ^ Projects the polynomial in the given dimension.

> instance Ord a => BsClass (Bs a) where
>   sharp k (B b) | sharpTreeL k b = Just L
>   sharp k (B b) | sharpTreeR k b = Just R
>   sharp _ _              = Nothing
>   proj k h (B b) = B $ projTree k h b

> instance Ord a => BsClass (Scaled a) where
>   sharp k (Scaled b _) = sharp k b
>   proj k h (Scaled b s) = Scaled (proj k h b) s
